function plot_heart_splines_combined()
    % Define the parameter range
    t = linspace(0, 2 * pi, 1000);
    [x, y] = heart_shape(t);

    % Calculate cumulative chordal length
    cumulative_lengths = cumulative_chordal_length(x, y);

    % Degrees for the splines
    degrees = [10, 40, 160];

    % Create a figure
    figure;

    % Original heart shape for reference
    plot(x, y, 'k--', 'LineWidth', 1);
    hold on;

    % Colors for each spline plot
    colors = ['r', 'g', 'b'];

    % Plot for different degrees
    for i = 1:length(degrees)
        n = degrees(i);
        
        % Parameterize t based on cumulative chordal length
        tck = linspace(0, cumulative_lengths(end), n);
        
        % Interpolate using cubic spline
        spline_x = spline(cumulative_lengths, x, tck);
        spline_y = spline(cumulative_lengths, y, tck);
        
        % Plotting the spline
        plot(spline_x, spline_y, 'Color', colors(i), 'LineWidth', 1.5);
    end
    
    % Formatting the plot
    title('Heart shape approximation with B-splines');
    legend({'Original', 'n=10', 'n=40', 'n=160'}, 'Location', 'Best');
    axis equal;
    hold off;
end

function [x, y] = heart_shape(t)
    x = 16 * sin(t).^3;
    y = 13 * cos(t) - 5 * cos(2*t) - 2 * cos(3*t) - cos(4*t);
end

function lengths = cumulative_chordal_length(x, y)
    dx = diff(x);
    dy = diff(y);
    lengths = [0, cumsum(sqrt(dx.^2 + dy.^2))];
end
